---
title: "R17.Migration"
output: pdf_document
---

# Load Packages
```{r}
library(tidyverse)
library(estimatr)
```

# Load Data
```{r}
load("Data/data_final.Rdata")
```

# Run Models
```{r}

output <- list()
for (g in 0:2){ # Subset for (3) groups (land reform intensity)
 
output[[g+1]] <- data.frame(covariates = "PRes_Change_7080", group = g, coef = NA, se = NA, n = NA)

  model <- lm_robust(as.formula(paste("PRes_Change_7080", paste(c("tx_Hi", "as.factor(DistN)"), collapse=" + "), sep = " ~ ")), 
                          data = data_final %>% filter(exp_cat2 == g), weight = weights)
  output[[g+1]][1,"coef"]<-model$coefficients["tx_Hi"]
  output[[g+1]][1,"se"]<-model$std.error["tx_Hi"]
  output[[g+1]][1,"n"]<-length(model$fitted.values)
  
}

# Plot Results

output_list <- do.call(rbind, output)  %>% 
  mutate(covnames = ifelse(covariates == "PRes_Change_7080", "% Change in Population Present, 1970-1980", NA))
output_list$covnames <- factor(output_list$covnames, 
                               levels = c("% Change in Population Present, 1970-1980"))       

ggplot(output_list, aes(x = reorder(covnames, desc(covnames)), col = as.factor(group), shape = as.factor(group), group = as.factor(group))) +
  geom_hline(aes(yintercept = 0), col = "red") +
  geom_point(aes(y = coef), position = position_dodge(width = 1)) + 
  geom_linerange(aes(ymin=coef - (qt(0.975, n)*se), ymax=coef + (qt(0.975, n)*se)), position = position_dodge(width = 1), size = 0.5) +
  geom_linerange(aes(ymin=coef - (qt(0.95, n)*se), ymax=coef + (qt(0.95, n)*se)), position = position_dodge(width = 1), size = 0.75) +
  coord_flip() +
  scale_color_manual(name = "Land Reform Intensity", labels = c("No LR", "LR < 50%, or adjacent", "LR >= 50%"), 
                     values = c("#252525", "#636363", "#969696"))+
  scale_shape_manual(name = "Land Reform Intensity", labels = c("No LR", "LR < 50%, or adjacent", "LR >= 50%"), 
                       values = c(16, 17, 15)) +
  theme_light() +
  theme(axis.text.x = element_text(angle = 90))+
  ylab("Treatment - Control Difference") + xlab("% Change in Population Present, 1970-1980") +
  theme(axis.text.y=element_blank())

ggsave(plot = last_plot(),
             "Figures/Afg6b_RobustMig.pdf",
             width = 7, height = 5)

```
